function [x, a] = SolveEquationSystemGaussCustom( A, b)
%Gets the solution of the equation system using gaussian elimination
C = sym([A b]);
N = size(A, 1);
M = size(A, 2);
for i = 1 : N - 1
    maxIndex = MaxElementIndex(A, i);
    C = ShiftRows(C, i, maxIndex(1));
    for k = i + 1 : N
        C = SubtractRows(C, i, k);
    end
end
result = sym(zeros(M, 1));
m = M;
while m >= 1
    sum = 0;
    p = M;
    index = N - (M-m);
    while p > m
       sum = sum + result(p) * C(m, p);
       p = p - 1; 
    end  
    if C(index, m) == 0
        result(m) = sym(strcat('C', num2str(m)));
    else
        result(m) = (C(index,M + 1) - sum) / C(index, m);
    end
    m = m - 1;
end

x = result;
a = C(:, 1 : M);

end

